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Abstract 

QCDFT is a multiscale modeling approach that can simulate multi-million 
atoms effectively via density functional theory (DFT). The method is based 
on the framework of quasicontinuum (QC) approach with DFT as its sole 
energetics formulation. The local QC energy is calculated by DFT with 
Cauchy-Born hypothesis and the nonlocal QC energy is determined by a 
self-consistent embedding approach, which couples nonlocal QC atoms to the 
vertices of the finite-elements at the local QC region. The QCDFT method 
is applied to a nanoindentation study of an Al thin film in the presence and 
absence of Mg impurities. The results show that the randomly distributed 
Mg impurities can significantly increase the ideal and yield strength of the 
Al thin film. 
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1. Introduction 

The ability to perform quantum mechanical simulations of materials prop- 
erties over length scales that are relevant to experiments represents a grand 
challenge in computational materials science. If one could treat multi-millions 
or billions of electrons effectively at micron scales, such first-principle quan- 
tum simulations could revolutionize materials research and pave the way to 
the computational design of advanced materials. 

In this paper, we propose a multiscale approach that is based entirely on 
density functional theory (DFT) and allows quantum simulations at micron 
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scale and beyond. The method, termed QCDFT tjPeng et all 120081 ). combines 
the coarse graining idea of the quasicontinuum (QC) approach and the cou- 
pling strategy of the quantum mechanics/molecular mechanics (QM/MM) 
method, and represents a major advance in the quantum simulation of ma- 
terials properties. It should be stated at the outset that QCDFT is not 
a brute-force electronic structure method, but rather a multiscale approach 
that can treat large systems - effectively up to billions of electrons. Therefore, 
some of the electronic degrees of freedom are reduced to continuum degrees 
of freedom in QCDFT. On the other hand, although QCDFT utilizes the 
idea of QM/MM coupling, it does not involve any classical/empirical poten- 
tials (or force fields) in the formulation - the energy calculation of QCDFT 
is entirely based on DFT. This is an important feature and advantage of 
QCDFT, which qualifies it as a bona fide quantum mechanical simulation 
method. 

Since QCDFT is formulated within the framework of the QC method, 
we shall give a brief introduction to QC in Sec. 2.1 to set up the stage of 
QCDFT. In Sec. 2.2, we briefly explain the local QC calculations. In Sec. 2.3, 
we introduce a DFT-based QM/MM approach that can treat the nonlocal 
QC region accurately and efficiently In Sec. 3, we apply QCDFT to the 
study of nanoindentation of an Al thin film in the presence and absence of 
Mg impurities. We present the nanoindentation results in Sec. 4 and finally 
our conclusions in Sec. 5. 



2. QCDFT Methodology 

2.1. Quasicontinuum Method 

The goal of the QC method is to mode l an atomistic system without ex- 



plicit ly treating every atom in the problem (ITadmor et allll996l ; IShenoy et al 



19991 ) . This is achieved by replacing the full set of N atoms with a small sub- 
set of N r "representative atoms" or repatoms (N r <C N) that approximate 
the total energy through appropriate weighting. The energies of individual 
repatoms are computed in two different ways depending on the deformation 
in their immediate vicinity. Atoms experiencing large variations in the de- 
formation gradient on an atomic scale are computed in the same way as in 
a standard atomistic method. In QC these atoms are called nonlocal atoms. 
In contrast, the energy of atoms experiencing a smooth deformation field on 
the atomic scale is computed based on the deformation gradient G in their 
vicinity as befitting a continuum model. These atoms are called local atoms 
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because their energy is based only on the deformation gradient at the point 
where it is computed. In a classical system where the energy is calculated 
based on classical/empirical interatomic potentials, the total energy .Etot can 
be written 

jV nl N loc 

£ff = X;^(R)+X;n i Ef*(G). (1) 

The total energy has been divided into two parts: an atomistic region of iV nl 
nonlocal atoms and a continuum region of N loc local atoms (iV nl +iV loc = N r ). 
The calculation in the nonlocal region is identical to that in atomistic meth- 
ods with the energy of the atom depending on the coordinates R of the 
surrounding repatoms. Rather than depending on the positions of neigh- 
boring atoms, the energy of a local repatom depends on the deformation 
gradients G characterizing the finite strain around its position. 

2.2. Local QC calculation with DFT 

In the local QC region, a finite element mesh is constructed with each 
repatom on the vertices of surrounding finite elements. The energy and force 
of each local repatom can be obtained from the strain energy density and 
the stress tensor of the finite elements that share the same repatom. More 
specifically, according to the Cauchy-Born rule, the deformation gradient G 
is the same for a given finite element, therefore the local energy density e and 
the stress tensor for each finite element can be calculated as a perfect infinite 
crystal undergoing a uniform deformation specified by G. In other words, one 
could perform a DFT-based energy/stress calculation for an infinite crystal 
by using periodic boundary conditions with the primitive lattice vectors of 
the deformed crystal, h, given by 

h^GHi, i = 1,2,3. (2) 

Here Hj are the primitive lattice vectors of the perfect undeformed crystal 
and Qo is the volume of the primitive unit cell. The deformed crystal can be 
derived from the perfect crystal via the deformation gradient G as shown in 
Fig.H) 

For the deformation gradient Gk associated with the fcth element, a pe- 
riodic DFT calculation can be performed to determine the strain energy per 
unit cell £' DFT (Gk). The Cauchy stress tensor can be defined as follows: 



^=oL In, ( 3 ) 



n ^ dh, 
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with being the volume of the deformed unit cell and hij denoting the 
component of the deformed lattice vector hj in Cartesian coordinate i. 

Once the strain energy E DFT (Gk) is determined, the energy contribution 
of the jth local repatom is given as 



M, 



Ef c ({G}) = J2^kE DFT (G k ), (4) 

k=l 

where Mj is the total number of finite elements shared by the jth repatom, 
and Wjk is the weight associated with the kth finite element for the jth local 
repatom. The force on the jth local repatom is defined as the gradient of the 
total energy with respect to its coordinate Rj OC - In practice, the nodal force 
on each finite element is calculated from the stress tensor of the finite element 



by using the Principle of Virtual Work (IZienkiewicz and Taylorl . 120001 ) . The 
force on the repatom is then obtained by summing the nodal force contribu- 
tions from each surrounding finite elements. 

2.3. Nonlocal QC calculation with DFT 

The nonlocal QC is modeled at the atomistic level with a QM/MM ap- 
proach. In a typical QM/MM calculation, the system is partitioned into two 
domains: a QM region and an MM region. In QCDFT, the QM atoms refer 
to the nonlocal repatoms and the MM atoms refer to the buffer atoms which 
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are the combination of both dummy atoms and local repatoms in QC termi- 
nology. The so-called dummy atoms are in the local region and are not inde- 
pendent degrees of freedom, but rather slaves to the local repatoms. In other 
words, the position of a dummy atom is determined by the finite element 
interpolation f rom t he relevant local repatom positions fjTadmor et all 11996 ; 



Shenoy et all Il999l ). The dummy atoms provide the appropriate bound- 



ary conditions for nonlocal DFT calculation while the energy of the dummy 
atoms is still treated with the Cauchy- Born rule, consistent with their status. 
The self-consistent embedding theory (jCholy et alll2005l ; IZhang and Lul . 12007 
Zhang et all 120081 ) is employed for the QM/MM calculations. More specif- 
ically, both the energy of the nonlocal atoms and the interaction energy 
between the nonlocal atoms and the buffer atoms are calculated by DFT. 
To simply the notation, we denote the nonlocal region as region I, and the 
buffer region as region II, as shown in Fig. [2j Typically, the buffer region 
consists of several atomic layers surrounding the nonlocal region. We asso- 
ciate each buffer atom in region II with a valence electron density (p at ) and 
a pseudopotential; both of them are constructed a p riori and remain fixed 
during a QM/MM simulation. (jZhang and Lul . 120071 ) The nonlocal energy 
E nl as defined in Eq. (1) can be expressed as 



E 



ni 



mm p i{E 



DFT 



[p 1 ; R 1 ] 



+ ^[p I ,p II ;R I ,R 11 ]}. 



(5) 



Here R 1 and R n denote atomic coordinates in region I and II respectively. 
The charge density of region I, p 1 , is the degree of freedom and is determined 
self-consistently by minimizing the nonlocal energy functional. The charge 
density of region II, p 11 , is defined as the superposition of atomic-centered 



charge densities p at via p (r) 



n int • 



Rj), which only changes upon 



the relaxation of region II ions. Eq£ is the interactio n energy between regions 
I and II computed by orbital-free DFT (OFDFT) flWang and Carter! . l2000l : 
Wang and Teterl Il992h . 

A basic ansatz of the nonlocal energy functional (Eq. (5)) is that p 1 must 
be confined within a finite volume (Q l ) that is necessarily larger than region 
I but much smaller than the entire QM/MM region. In addition, since some 
terms in the for mulation of Eq. (5) c ould be more efficiently computed in 



reciprocal space (jZhang and Lul . 120071 ). we also introduce a volume Q B over 



which the periodic boundary conditions are applied. The periodic box fl B 
should be large enough to avoid the coupl ing errors induced by the imple- 
mentation of periodic boundary condition (jZhang and Lul . 120071 ). 
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Figure 2: The schematic diagram of domain partition in QCDFT with a dislocation in 
Al lattice as an example. The black and white spheres represent the nonlocal and buffer 
atoms, respectively. The dotted box represents J7 1 and the solid box represents the periodic 
box f2 B . The volume fi 1 and f2 B is 2.8 A and 8 A beyond the region I in ±x and ±y 
directions, respectively. 



The interaction energy, Eq^, formulated by OFDFT is defined as follow- 



ing: 



E^[ P \ p n ; R\ R 11 ] = £ OF [p tot ; R tot ] - ^orfp 1 ; R 1 ] 

— E OF [p u ; R n ], 



(6) 



where R tot = R 1 |J R 11 and p tot = p 1 + p n . In addition to its computational 
efficiency, OFDFT allows Eq. (9) to be evaluated over Q l rat her than over 
the entire QM/MM syste m as Eq. (9) appears to suggest (jCholy et al. 



20051 ; IZhang and Lul . 120071 ). This significant computational saving is due to 
the cancellation in evaluating the first and second term of Eq. (9), and it 
is rendered by the orbital-free nature of OFDFT and the localization of p 1 . 
A single-particle embedding potential p em b(r) can be defined as a functional 
derivative of the interaction energy with respect to p 1 



Pemb( r ) = 



5p l 



(7) 



which represents the e ffective potential that region I elec t rons f eel due to the 



presence of region II; (jCholy et all 120051 ; IZhang and Lul . 120071 ) it is through 
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Figure 3: The solid charge density of the perfect Al lattice along [100], [110] and [111] 
directions obtained by the periodic DFT calculation and the superposition of the fitted 



Pemb(r) that the QM/MM coupling is achieved quantum mechanically at 
the level of OFDFT. The embedding potential provides rigorous boundary 
conditions for p 1 and is updated self-consistently during the minimization of 
the nonlocal energy functional. 

Since p 11 is a key quantity for the accurate calculation of the interaction 
energy and the embedding potential, it is crucial to construct an appropriate 
representation of p n . In fact, the construction of an appropriate charge 
density distribution in region II represents a common challenge to many 



QM/MM methods. (ILin and Truhlarl . 120071 ) In this paper, we represent p n 



as a superposition of spherical atomic-like charge densities centered on each 
ions in region II, which is a good approximation for metallic systems. Ideally, 
the constructed p n (r) should reproduce the bulk (or solid) charge density 
obtained by a DFT calculation of the perfect lattice. That is to say, one 
needs to determine p at (r) by minimizing the function 

/ [p n (r) - p soUd (r)] 2 dr (8) 
•/v u 

with p n (r) = p at (r — R^)- Here V u represents the volume of the unit cell, 
and p solld is the solid charge density obtained by a periodic DFT calculation 
for the perfect reference system. The summation of p includes all the ions 
which have contribution to the charge density in the unit cell. 

In this paper, we employ the param eterized multiple Slater-type orbitals 



(MSTO) flClementi and Roettil . 11974 ) for the expansion of p at (r). With 



MSTO, the atomic wave function $ of many-electronis is the superposition 
of all relevant atomic orbitals: $(r, 9,<p) = ^ Q0j(r, 9, <p), where q is the 



7 



weight of orbital i in the expansion and the ith atomic orbital can be written 

as 



Ur,e^) = Ar n - 1 e^ r Y l m {9^) ) (9) 

where n, I, m are the principal, angular momentum and magnetic quantum 
number of the orbital. Yi m (9, ip) is spherical harmonic function and ( is 
related to the effective charge of the ion. A is a normalization constant and 
is expressed as A = (2£)™ +1 / 2 / a/ (2n)\. With this expansion, the atomic- 
centered charge density can be calculated as 

/7T /*27T 
/ \<S>*(r,9,<p)$(r,9,<p)\d9d<p. (10) 

The parameters Cj and ( are determined by minimizing Eq. ([8]) with the 
constraint of preserving the correct number of valence electrons. In Fig. [31 
we present the solid charge density p solld determined from p at and from the 
periodic DFT calculation for a perfect Al lattice. It can be seen that the 
constructed p solld reproduces very well the solid charge density calculated by 
DFT calculations for the same perfect lattice. 



3. Computational details 



3.1. Model setup 

The present QCDFT approach is applied to nanoindentation of an Al 
thin film resting on a rigid substrate with a rigid knife-like indenter. The 
QCDFT method is appropriate for the problem because it allows the mod- 
eling of system dimensions on the order of microns and thus minimizes the 
possibility of contaminating the results by the boundary conditions arising 
from small model sizes typically used in MD simulations. The reason we 
chose this particular system is because the re exists a good kine t ic en ergy 
functional and an excellent EAM potential (jErcolessi and Adamd . Il994j) for 
Al. In this paper, we have resc aled the "force-matching" EAM potential of 
Al ( jErcolessi and Adamd . If 994l ) so that it matc hes precisely the DFT value 



of the lattice constant and bulk modulus of Al (IPeng et all 120081 ) 



The crystallographic orientation of the system is displayed in Fig. HI The 
size of the entire system is 2 /im x 1 /iin x 4.9385 A along the [111] (x 
direction), the [110] (y direction), and the [112] (z direction), respectively. 
The system is periodic in the z-dimension and has the Dirichlet boundary 
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Figure 4: Schematic representation of the nanoindentation of Al thin film showing the 
relevant dimensions and orientations. 



conditions in the other two directions. The entire system contains over 60 
million Al atoms - a size that is well beyond the reach of any full-blown 
quantum mechanical calculation. The unloaded system is a perfect single 
crystal similar to the experimental situation. The film is oriented so that the 
preferred slip system (110) {111} is parallel to the indentation direction to 
facilitate dislocation nucleation. The indenter is a rigid flat punch of width 
25 A. We assume the perfect-stick boundary condition for the indenter so 
that the Al atoms in contact with it are not allowed to slip. The knife-like 
geometry of the indenter is dictated by the pseudo two-dimensional (2D) na- 
ture of the QC model adopted. Three-dimensiona l QC models do exist and 



can b e implemented in QCDFT (IKnap and Ortizl . 120031 ; lHayes et all 120051 . 



20061 ). We chose to work with the pseudo-2D model in this example for its 



simplicity. The prefix pseudo is meant to emphasize that although the anal- 
ysis is carried out in a 2D coordinate system, the out-of-plane displacements 
are allowed and all atomistic calculations are three-dimensional. Within this 
setting only dislocations with line directions perpendicular to the xy plane 
can be nucleated. 

The simulation is performed quasistatically with a displacement control 
where the indentation depth (d) is increased by 0.1 A at each loading step. 
Because D FT calculations are m uch more expensive than EAM, we use EAM- 
based QC flTadmor et all Il999bh to relax the system for most of the loading 



9 



steps first. At d = 2.0,3.0,5.5,6.0,7.0,7.1,7.5 A, the corresponding EAM 
configurations are further relaxed by QCDFT. The QCDFT loadings are car- 
ried out after d = 7.5 A starting from the full relaxed EAM-QC configuration 
of previous loading step, until the onset of the plasticity occurs at d — 8.2 
A. Such a simulation strategy is justified based on two considerations: (1) 
an earlier nanoindentation study of the same Al surface found that the onset 
of plasticity occurred at a smaller load with EAM-base d local QC calcula - 
tions comparing to OFDFT-based local QC calculations (IHayes et al.L 120051 ) . 
The result was obtained by a local elastic stability analysis with EAM and 
OFDFT calculations of energetics and stress. This suggests that we will not 
miss the onset of plasticity with the present loading procedure by performing 
EAM-QC relaxations preceding QCDFT. (2) Before the onset of plasticity, 
the load-displacement response is essentially linear with the slope determined 
by the elastic properties of the material. In other words, only two QCDFT 
calculations are required to obtain the linear part of the loading curve. 



Figure 5: Schematic diagram of the randomly distributed Mg impurities in the Al thin film. 
The red spheres and blue pentagons represent nonlocal Al and Mg atoms, respectively. 
The green triangle represents Al buffer atoms. The dimensions are given in A . 

We also study the effect of Mg impurities on the ideal strength and incip- 
ient plasticity of the Al thin film. In the calculations, five Mg impurities are 
introduced randomly below the indenter, as schematically shown in Fig. [5j 
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The results of the randomly distributed Mg impurities are referred as ran- 
dom, distinguishing from the results of the pure system, referred as pure. 
At d = 3.0, 6.0, 7.5 A, the random results are obtained after full relaxations 
of the pure Al system. The QCDFT loading is carried out after d = 7.5 A 
starting from the full relaxed configuration of a previous loading step, until 
the onset of the plasticity occurs at d — 8.1 A. 

The p arameters of t he OF DFT density-dependent kernel are chosen from 
reference (jWang et all Il999l ). and Al ions are represented by the Goodwin- 



Needs-Heine local pseudopotential ( Goodwin et al. . 199Cll ). The high kinetic 
energy cutoff for the plane wave basis of 1600 eV is used to ensure the con- 
vergence of the charge density. For the nonlocal calculation, the grid density 
for the volume Q 1 is 5 gridpoints per A. The n 1 box goes beyond the nonlocal 
region by 8 A in ±x and ±y directions so that p l decays to zero at the bound- 
ary of Q l , as shown in Fig. [2j The relaxation of all repatoms is performed 
by a conjugate gradient method until the maximum force on any repatom is 
less than 0.03 eV/A. 



4. Results and Analysis 

The load-displacement curve is the typical observable for nanoindenta- 
tion, and is widely used in both experiment and theory, often serving as a link 
between the two. In particular, it is conventional to identify the onset of in- 
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present work, the load is given in N/m, normalized by the length of the 
indenter in the out-of-plane direction. 

For pure Al, the load-displacement (P — d) curve shows a linear relation 
followed by a drop at d — 8.2 A, shown by the dashed line in Fig. [6J The 
drop corresponds to the homogeneous nucleation of dislocations beneath the 
indenter - the onset of plasticity. A pair of straight edge dislocations are 
nucleated at x=±13 A, and y=-50 A. In Fig. [3, we present the out-of-plane 
(or screw) displacement u z of the nonlocal repatoms. The non-zero screw 
displacement of the edge dislocations suggests that each dislocation is disso- 
ciated into two 1/6 <112> Shockley partials bound by a stacking fault with 
a width of about 19 A. The activated slip planes are those {111} planes that 
are adjacent to the edges of the indenter. The slope for the linear part of 
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the curve is 27.1 GPa, which is less than the shear modulus /i=33.0 GPa 
and C44 = 29.8 GPa. The critical load, P cr for the homogeneous dislocation 
nucleation is 18.4 N/m, corresponding to a hardness of 7.3 GPa (the critical 
load normalized by the area of the indenter), which is 0.22 \x . The drop in 
applied load due to the nucleation of dislocations is AP = 6.8 N/m , agree - 



ing with the load drop estimated by the elastic mode hlTadmor et al 
which is AP = 7.7 N/m. 
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Figure 6: Load-displacement plot for the nanoindentation of the Al thin film with a rigid 
rectangular indenter: pure Al (red squares) and randomly distributed Mg impurity system 
(green circles). The corresponding lines are the best fit to the data points. 

For randomly distributed impurities in the Al thin film, the load-displacement 
curve shows a linear relation up to a depth of 8.0 A, followed by a drop at 
d = 8.1 A, as shown by the solid line in Fig. [6j The slope of initial linear part 
of the load-displacement curve is 26.7 GPa, rather close to the correspond- 
ing pure Al value. The maximum load in linear region is P^P" = 19.2 N/m, 
corresponding to a hardness of 7.6 GPa, which is 0.3 GPa greater than the 
pure Al system. A pair of Shockley partial dislocations is nucleated at x=-13 
A, y=-25 A and x=13 A, y=-22 A respectively as shown in the right panel 
of Fig. The drop in the applied load due to the dislocation nucleation is 
5.9 N/m. The estimated load drop by the elastic model is AP = 7.6 N/m. 
The smaller drop of the load for the random case than the elastic model is 
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probably due to the pres ence of the Mr impur ities, which is not accounted 



for in the elastic model (ITadmor et all Il999bl ). The fact that the critical 



load and the hardness of the Al-Mg alloy are greater than that of the pure 
Al system demonstrates that the Mg impurities are responsible for the solid 
solution strengthening of the Al thin film. The presence of Mg impurities also 
hinders the formation of full edge dislocations and as a result, only partial 
dislocations are nucleated and they are pinned near the surface as shown in 

Fig.m 




Figure 7: The out-of-plane displacement u z obtained from the pure (left) and with Mg 
impurities (right) QCDFT calculations. The circles represent the repatoms and the dis- 
placement ranges from -0.4 (blue) to 0.4 (red) A. 



Finally we point out the possibility that the emitted dislocations may 
be somewhat constrained by the local/nonlocal interface from going further 
into the bulk. Because the critical stress to move an edge dislocation in Al 
is vanishingly small ( 10 -5 /i) comparing to that to nucleate a dislocation 
(lCT 1 ^), a small numerical error in stress could easily lead to a large dif- 
ference in the equilibrium dislocation position. The four-order-of-magnitude 
disparity poses a significant challenge to all atomistic simulations in predict- 
ing dislocation nucleation site, QCDFT method included. One can only hope 
to obtain a reliable critical load for the incipient plasticity, rather than for 
the equilibrium position of dis locations. The same problem has been ob- 
served and discussed by others (ITadmor et all Il999al ). However, despite the 



problem, the dramatic difference observed in the two panels of Fig. 7 unam- 
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biguously demonstrates the strengthening effect of Mg impurities. Therefore 
the conclusion is still valid. 

5. CONCLUSION 

In summary, we propose a concurrent multiscale method that makes it 
possible to simulate multi-million atoms based on the density functional the- 
ory. The method - QCDFT - is formulated within the framework of the 
QC method, with DFT as its sole energy input. The full-blown DFT and 
DFT-based elasticity theory would be the two limiting cases corresponding 
to a fully nonlocal or a fully local version of QCDFT. The QCDFT method 
is applied to nanoindentation of an Al thin film in the presence and absence 
of randomly distributed Mg impurities. The Mg impurities are found to 
strengthen the hardness of Al and hinder the dislocation nucleation. The 
results suggest that QCDFT is a promising method for quantum simulation 
of materials properties at length scales relevant to experiments. 
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